Precession— tracking coordinates for simulations of compact— object— binaries 
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Binary black hole simulations with black hole excision using spectral methods require a coordi- 
nate transformation into a co-rotating coordinate system where the black holes are essentially at 
rest. This paper presents and discusses two coordinate transformations that are applicable to pre- 
cessing binary systems, one based on Euler angles, the other on quaternions. Both approaches are 
found to work well for binaries with moderate precession, i.e. for cases where the orientation of the 
orbital plane changes by <C 90° . For strong precession, performance of the Euler-angle parameteri- 
zation deteriorates, eventually failing for a 90° change in orientation because of singularities in the 
parameterization ("gimbal lock"). In contrast, the quaternion representation is invariant under an 
overall rotation, and handles any orientation of the orbital plane as well as the Euler-angle technique 
handles non-precessing binaries. 
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I. INTRODUCTION 

Gravitational waves offer an exciting new observational 
window into the universe. With the second generation of 
gravitational wave detectors such as Advanced LIGO and 
Advanced Virgo commencing observations in 2015 [l|, it 
is extremely important to develop a detailed picture of 
the gravitational physics of the most likely sources. A 
very promising source of gravitational waves are inspiral- 
ing and merging binary black holes Q- Because of the 
weakness of the gravitational wave signal, matched filter- 
ingis necessary to pick out the waveform from the noise 
[sl, |j|. Constructing such templates, in turn, requires di- 
rect numerical integration of Einstein's equations for the 
late inspiral, merger and ringdown phase of the coalesc- 
ing compact object binary; see, e.g. [5[. Since 2005, start- 
ing with the seminal work of Frans Prctorius [y] , many 
groups have successfully simulated binary black hole sys- 
tems using a variety of different techniques. For recent 
overviews of the state of the field, see 0, Q- 

Compact object inspirals fall into two categories: non- 
precessing and precessing. While the non-precessing 
aligned-spin systems arguably represent an important 
subspace of all binary black hole systems, the more gen- 
eral case features arbitrary spin orientations. In this 
non-symmetric situation, the interaction of the orbital 
angular momentum and the black holes' spins leads to 
precession of the orbital plane, changing its orientation 
by as much as 180 degrees. 

Precession modulates the gravitational waveform. 
Therefore, it is crucial to explore these strongly pre- 
cessing systems. Furthermore, precessing systems allow 
the study of gravitational dynamics in an underexplored 
regime, providing a new opportunity for comparing nu- 
merical relativity to various analytic approximations like 
post-Newtonian (see e.g. [3, [i3|) and effective-one-body 
theory (see e.g. |lll4l4| ). The numerical simulations can 



be used both to test the accuracy of the analytic treat- 
ments and to calibrate them, in some cases, thus improv- 
ing their accuracy |15l4l7l | . Furthermore, one can attempt 
to reproduce numerically predictions from analytic com- 
putations such as transitional precession [ISj . which is 
known from PN theory but has not yet been observed in 
numerical simulations. 

Numerical simulations of precessing binary black holes 
have already been undertaken; for example |19l - [26| . 
Given the vastness of parameter space and the need for 
simulations lasting at least 10 orbits - possibly lOO's of 
orbits - to optimally exploit gravitational wave detec- 
tors |5l. l27H3l| . a lot of extra work remains to be done. 

The Spectral Einstein Code SpEC [33 is allows effi- 
cient and accurate simulations of binary black holes; see 
e.g. [3l|, I33l-l39|. This code applies black hole excision 
and uses time-dependent coordinate mappings to rotate 
and deform the computational grid such that the ex- 
cision regions remain inside the black hole horizons at 
all times. For non-precessing inspiralling binaries, these 
coordinate mappings are described in detail in previous 
work [H, 111,13 • 

The purpose of the present paper is to develop coordi- 
nate mappings that are able to follow a precessing com- 
pact object binary through the inspiral, even for strongly 
precessing systems. Wc present two different approaches. 
The first one is based on Euler angles; it works well for 
moderate precession, but fails when the orientation of the 
orbital plane changes by 90 degrees or more. The second 
approach is designed to avoid the deficiencies of the Euler 
angle parameterization. By using quaternions, we devise 
coordinate mappings that work for any change of orienta- 
tion of the orbital plane with a performance comparable 
to the earlier non-precessing techniques. The techniques 
developed here have already been used in [4l|, [d^l 

This paper is organized as follows. Section |lT] de- 
scribes the computational setup of SpEC in more detail 



fScc. lII A| . and develops the coordinate mappings based 
on Euler angles (Sec. IIIBp and quaternions (Sec. IIICp . 
Section IIIII presents a sequence of numerical results ob- 
tained with both approaches, starting from Newtonian 
and post-Newtonian test-cases to simulations of binary 
black holes with numerical relativity (NR). We summa- 
rize our results in Sec. IIVI 



II. METHODS AND TECHNIQUES 

A. Dual frames and control systems 

As described in Scheel et al. [4j], SpEC utilizes a 
dual-frame approach to simulate compact object bina- 
ries. Einstein's equations are written down in an asymp- 
totically non-rotating coordinate-system x"" = (t, a;*), re- 
ferred to as the "inertial frame" , and all tensors are rep- 
resented in the coordinate basis of this frame. In the 
inertial frame, tensor components remain finite even at 
large separation. The computational grid is specified in 
"grid coordinates" x" = (t,a;'). The collocation points of 
the spectral expansion are at constant grid coordinates, 
and numerical derivatives are computed with respect to 
these coordinates. The two coordinate frames share the 
same time-coordinate 



t = t. 



(1) 



The spatial coordinates of the two frames are related by 
a coordinate transformation 



;A^W) 



(2) 



which depends on a set of parameters X^{t) to be dis- 
cussed in detail later. The coordinate transformation 
Eq. ([2|) maps the grid-coordinates into the inertial frame 
such that the excision surfaces (coordinate spheres in the 
grid-frame) are mapped to locations somewhat inside the 
apparent horizons of the black holes in the inertial frame. 
In the original work [4^, this coordinate transforma- 
tion was taken as the composition of a rotation about the 
z-axis and an overall scaling of the coordinatefu, 
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(3) 



In this simple case, the map A^(i) = {a{t),tp(,t)} de- 
pends on two parameters: the scale factor a{t) and the 
rotation angle ip{t). The map parameters A^(i) are cho- 
sen dynamically during the simulation, such that the map 
tracks the actual motion of the black holes. This can be 
accomplished by introducing a set of control-parameters 
Q^, such that 



1. Q'^ = if the mapped excision spheres are at the 
desired location in inertial coordinates. 

2. Under small variations of the mapping parameters 
around their current values, the control-errors sat- 
isfy 



dQ^ 



dX'' 



'6it 



(4) 



Af =Af'(t) 



While not strictly required, Eq. Q allows one to write 
down uncoupled feedback control equations for the A'^(i). 
In the special case of a linear, uncoupled system, this 
reduces to Qf^ = X^^rget ~ ^'^■ 

For black holes orbiting in the xy-plane, Eq. (jH]) suf- 
fices to keep the excision boundaries inside the inspiraling 
black holes, resulting in successful simulations of inspiral- 
ing BH-BH binaries in Ref. [4J|. Subsequently, the map 
was refined to avoid a rapid inward motion of the outer 
boundary [3J], to adjust the shapes of the mapped exci- 
sion boundaries to more closely conform to the distorted 
apparent horizons |37l . l38l . |45|, and was generalized to 
unequal mass binaries |36| . Hemberger et al [40| summa- 
rizes these maps, and introduces further mappings that 
are needed during the merger phase of the black hole 
binary. 

The purpose of the present paper is the development of 
coordinate mappings that can handle precessing binaries. 
Because in general the center of mass will move (e.g due 
to asymmetric GW emission) , these coordinate mappings 
must also allow for a translation of the binary. Rotation 
and translation couple to each other and must therefore 
be dealt with simultaneously. The questions addressed in 
this paper are therefore (1) determination of a suitable 
coordinate mapping for precessing, translating binaries, 
(2) suitable choice of mapping parameters A^, and (3) 
derivation of control-parameters Q^. Specifically, we will 
discuss below two generalizations of Eq. (|3]), one based 
on Euler-angles and one based on quaternions. We will 
show that the Euler-angle representation suffers from sin- 
gularities when the inclination of the orbital plane passes 
through 7r/2, and we will demonstrate that the quater- 
nion representation fixes these problems. 



B. Euler angle representation 

In the general case where the orbital plane precesses, 
we use a mapping that composes a scaling a(t), a rotation 
R(t) and a translation T{t). The mapping is given by 

x = a{t)R{t)x + f. (5) 

A rotation matrix can be specified by Euler angles. 



'cos COS ip — COS (j) sin ip + sin sin 9 cos ip sin </> sin V' + cos (p sin 6* cos ip 
R =^ I cos 6 sin -0 cos </> cos V" + sin sin sin ?/; — sin cos f/; + cos sin sin ip 
— sin 6* sin 6 cos 6* cos 6 cos 61 



(6) 



where (j) is the roll angle around the x-axs, 9 is the pitch 
angle around the y-axis, and ip is the yaw angle around 
the z axis and we have suppressed the explicit time- 
dependence. 

For our application, the desired locations of the black 
holes lie parallel to the x-axis, i.e. the black holes are 
at grid coordinates [0^,0^,0^) and (c|;,c^,c^). It is 
straightforward to show that for these two points a ro- 
tation about the x-axis is degenerate with a translation 
because only the location of the black holes is important. 
Therefore, we can set (j){t) = so thal[j 



cos 9 cos ip 

R = I cos 6 sin ip 

— sin 9 



' sin ip sin 9 cos ip^ 
cos ip sin 9 sin ip 
cos( 



(7) 



Thus the mapping in Eq. ([S|) will have six parameters in 
this case, a scaling a(t), a pitch angle (rotation about y- 
axis) 9(t), a yaw angle (rotation about z-axis) ip{t), and 
a translation {T^ {t),T^{t),T^{t)). 

The goal of the scaling-rotation-translation map is to 
keep the horizons of the two black holes centered on the 
excision surfaces. As the binary evolves, the map param- 
eters need to be adjusted by the control system. To derive 
the control parameters Q'^, c.f. EqQ, one can consider 
perturbations of the mapping parameters around their 
current values. Let A'' = {a,9,ip,f = (T^,T^,T^)} be 
the current imperfect mapping parameters at some time 
during the evolution. Furthermore, denote the desired 
parameters \^ = A'^ + 5X>' = {aQ,9o,ipo,Tf^ ,T^ ,Ti}- 
Finally, let xa and xb denote the current location of 
the center of black hole A and B, respcctivclj|j, and let 
CA and cb denote the desired location of the black hole 
centers; i.e., the centers of the excision spheres. For con- 
venience we also define the vectors X = xa — xb and 
C = ca — ca- The target mapping Aq is such that the 
points Ca, cb are mapped onto the inertial frame position 
of the black holes, ca = xa, cb = xb' 



XA.B = aoR{9o, 1po)cA,B + Tq. 



(8) 



^ Note that in Ref . |43|| , the equations give the transformation from 
the inertial coordinates to the grid coordinates, the rotation angle 
is <j) instead of ^, and the scale factor a is the inverse of the scale 
factor in this paper. 

^ Note that if the motion is confined to the x — y plane, the pitch 
will remain fixed at 9 = and we recover the rotation matrix in 
Eq. ©. 

^ The precise definition of "center" is not important; wc shall use 
the coordinate point around which the coordinate radius of the 
apparent horizon has vanishing / = 1 multipoles. 



r 



Rewriting this equation in terms of the current mapping 
and the grid location of each black hole (i.e xa.b) yields 

aR{9, ip)x,+f = ia+6a)R{9+69, i'+Si')c^+f+Sf. (9) 

where i — A, B. 

Equation ^ represents six equations for the six un- 
knowns 5\^ . Solving this system of equations to leading 
order in the perturbations yields 



6a = a\ — — 
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5ip = 
ST^ = -^ {6t^ cos 9cosip- 5t^ sin iP 



cos 9 X^ 
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ST^ = -^ {St^ cos sin iP + St^ cos ip 



6t^ 



where 
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6t-' 
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sin 9 sin ip) , 

-5t^ sin 9 + St^ cos ( 



Ca^b - c'bXa + cVXy + c'X' 



ot = caVb - CbVa - c^X 
- c^Xnan6l, 

(5t^ = c\zB - c%ZA -f cyxy\.&ni 
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(10c) 

(lOd) 

(lOe) 
(lOf) 

(11a) 
(lib) 
(lie) 



Furthermore, we have assumed that the centers of the 
excision surfaces are aligned parallel to the x-axis so that 

c'a — ^B— ^ ^"^^ c^ = c^ = c^. 

Perhaps surprisingly, the 8\^ given by Eq. (|10al llOfj) 
are the desired control parameters Q^ . This can be seen 
as follows. For a perfect map, A^ = Ag . i.e. 6\^ = Q^ = 
0. Moreover, by definition 



^ - A (A- - A-) - ^ - ^ - -5- (12) 



Thus (5A^ defined by Eq. (|10a| - [T0^ satisfy the conditions 
for Q^ outlined in section III Al 

The Euler angle prescription as described above is ad- 
equate for describing rotations that are close to the x — y 
plane, and has been used for the SpEC simulations pre- 
sented in [45|. However, the Euler angle prescription car- 
ries with it an inherent coordinate singularity that causes 
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FIG. 1. Typical behaviour of the Euler angles and their 
derivatives for a nearly polar orbit inclined at 85 degrees with 
respect to the x — y plane 



a breakdown of the control system for high inclination 
angles. Firstly, note that Eq. ()10c|) shows that 5-0 will 
diverge when 6 = ^, which would lead to the breakdown 
of a feedback control system. Further, since SpEC uses 
proportional-derivative-control or proportional-integral- 
derivative-control [401, we must examine the behaviour 
of the derivatives of the Euler angles, 9, ip. Notice that 
we can relate the angular velocity to the derivatives of 
the Euler angles simply by the relationship 



where A is the Euler angle rates matrix: 



' cos?/" cos — sin-0 0"* 
A = I sin tjj cos 9 cos ip 
-sin 6* 1, 



(13) 



(14) 



The black holes move on regular trajectories with a slowly 
varying orbital frequency Co in inertial coordinates; there- 
fore, the left-hand side of Eq. ([T^ is continuous. How- 
ever, |detA| = |cos(0)|, so that for ^ = § the time- 
derivatives of the Euler angles will diverge since, by 
Cramer's rule, the inverse oi A scales as 1/det^. 

There is another way to envision the divergence of the 
derivatives of the Euler angles. Consider the unit vector 
in the direction connecting the centers of the two compact 
objects in inertial coordinates, u = ,5"~5^, . Before, we 

■> ' \xb-xa\ ' 

considered "0, 9 as parameters in a mapping. Let us 
now consider them as spherical polar coordinates that 




cos tp cos t 
sin -0 cos ( 
sin( 



(15) 



We can immediately derive the expressions for 9, "0 as 
functions of Mr , w„, it,: 







cos 9' 
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cos 6* 



w? tan^ 



(16) 
(17) 



From these equations it is obvious that the derivatives of 
9 and behave abnormally when ii moves across one of 
the poles at uniform velocity. In fact, 9 does not exist, 
and its second derivative diverges (see the top panels of 
Fig. [T]). Meanwhile, diverges: letting 69 = ^ — 9 we 
can write, for (50 <C 1: 



ip oc 



U! 



(18) 



This behaviour is demonstrated clearly in Figure 
[TJ where the derivatives of both Euler angles angles 
demonstrate sharp and nearly discontinuous features. 

This is the fundamental reason why Euler angles are 
not a suitable parametrization of rotations: there ex- 
ist situations, which we would like to study, when their 
derivatives grow extremely fast numerically. 



C. Rotation— invariant Quaternion representation 

The origin of the break-down of the Euler angle rep- 
resentation lies in its reliance of a preferred coordinate 
system, which is implicit in the adoption of Euler angles. 
The physics of compact object inspirals is invariant under 
rotations of the spatial coordinates. Ideally, the numer- 
ical methods used to describe such a system should also 
be invariant, and should work equally well independent 
of the orbital plane of the black holes. 

The singularities in the Euler angle representation 
arise from a poor choice of representation of the rota- 
tion group, which relied on preferred directions in space 
(namely the coordinate axes). Therefore, a suitable rep- 
resentation must be independent of any special direc- 
tions. We employ quaternions to represent rotations and 
build up the overall rotation from a sequence of infinites- 
imal rotations. 

It should be noted that a similar construction can be 
done with a different paramtetrization of rotations. See 
the Appendix for an example using orthogonal infinites- 
imal rotation matrices. 



* 9 here means the angle to the xy plane rather then the angle to 
the z-axis. Therefore the following equations differ slightly from 
standard spherical polar coordinates. 



1. Quaternion algebra 

Quaternions are an extension of the complex numbers, 
with three imaginary units i,j, and k, obeying 



i^ = f = fc2 = ijk = -1, 



(19) 



as weU as certain further muhiphcation rules. A quater- 
nion q has the form 



q = qa + qii + q2J + q^k, go, ■ • • , 93 e 



(20) 



This is conveniently written as q = {qQ,q), where 
q~ (gij 92, ga)- Addition and scalar multiplication are 
defined in analogy with complex numbers. With the 
structure introduced so far, the set of quaternions 



implies that v' has indeed a vanishing real part. Equa- 
tion (|27p is equivalent to v' = RqV, with rotation matrix 

/ql + 2ql-q^ 2{qiq2~qoq3) 2{qaq2 + qiq3)\ 

Rq = \ 2{qiq2 + qoq3) ql + 2ql- q"^ 2(^293-9091) (28) 

\2(9i93-9o92) 2(go9i+9293) ql + '^ql- q^ ) 

We can now rewrite Eq. ([5]) in quaternion language, 
x = aqxq* +T. (29) 

Here, x,x,T have been promoted to quaternions; e.g., 
T= (0,f). 

Our next task is to derive equations that determine the 
time evolution of q as well as the control parameters Q'^. 



= {qo + qii + q2J + q^Mqi £ K} 



(21) 



is a 4-dimensional vector space over the real numbers. 
Multiplication is defined by 

gp= (Po9o -p-9,Po9 + 9oP + P X g), (22) 

where q ■ q and q x p are the standard Euclidean dot 
and cross products respectively. Complex conjugation is 
given by 

9* = (90, -9^. (23) 

It follows that the multiplicative inverse is given by 
.-1 _ 1* 



where the norm |q| satisfies 



\q\^ = qq*=ql + q^- 



(24) 



(25) 



Restricting our attention now to the set of all unit 
quaternions 5*15(1) = {q <E M, \q\ = 1}, it is easy to show 
that Sp{l) is isomorphic to SU{2) where SU{2) is the 
group of all 2 X 2 unitary matrices with unit determi- 
nant [43. SU{2) is a double cover of the rotation group 
50(3), which means that unit quaternions do represent 
rotations. 

Unit quaternions are related to rotations in the follow- 
ing manner. Let n be a unit-vector, and define 



o = cos — , n sm — 

^ ^ 2 2^ 



(26) 



for some angle 6. The quaternion q rotates a vector v 
into the vector z/', around the axis fi by angle 9 in the 
right-handed sense via 



V = qvq 



(27) 



In this equation, 3- vectors are to be promoted to quater- 
nions by the rule v = {0,v), v' — (0,w'),and |q| — 1 



2. Quaternion kinematics 

In this section we derive the differential equation 
obeyed by the rotation quaternion q. Consider a time- 
dependent unit-quaternion q(t) : M. -^ Sp{l). The 
derivative is defined by 



.^g(t + M-,(t) 

h^O h 



(30) 



We write the rotation q{t + h) at time t + h, as a product 
of q{t) and a quaternion u representing an infinitesimal 
rotation. 



q{t + h)=uq{t). 



(31) 



The quaternion u is easily obtained by expanding the 
right hand side of Eq. (P5| to first order in 6 using cos | « 
1, sin I w I: 



u^{l,n-) = l + dq/2 



(32) 



with Sq = {0,n6). If the rotational velocity is to in in- 
ertial coordinates, then n ~ uj and 69 = u h, so that 
Sq = u) h. Thus we can write q{t + h) = {I + 8q/2)q{t). 
Substituting, 

. ,. {I + 8ql2)q{t)-q{t) 1 

q = hm = -U3 q. (33 

h^a h 2 ^ 

Noting that in grid coordinates the angular velocity is 
O = q*ujq, we finally obtain |47lpl 



q = 2^^- 



(34) 



This reference uses the opposite convention of the one adapted 
here: Q is the angular velocity in the fixed frame whereas ui is 
the angular velocity in the rotating frame. 



3. Quaternion control system 

The goal of the sector of the control-system for rota- 
tions is to keep the vector X parallel to the vector C. 
The misalignment between them can be measured by the 
rotation needed to make these vectors parallel: 



Qr = 



C xX 



(35) 



where the subscript 'R' indicates that this quantity is 
of relevance for rotationsQ The control-system needs to 
adjust the angular velocity ft such that Qji pa 0. As long 
as the control-system works, this instantaneous rotation 
is small, and therefore, non-commutativity of rotations 
can be neglected. This suggests to control the angular 
velocity in the moving frame D, based on the control- 
parameter Qn. 

We proceed as follows: We measure Qr regularly dur- 
ing the BBH evolution, and compute its first and second 
time-derivatives. As in earlier work [43| (and in many 
papers since [33l - l36i ISa. l39l |44| ) . we use this to reset the 
third time-derivative of the mapping-parameters that de- 
termine the rotation. These parameters are the second 
time-derivative of fl{t); thus, we choose i7(i) such that 
it has constant second time-derivative. We periodically 
reset this constant using the equation 






aQu 



dQn (PQr 

l^—rr+l- 



dt 



dt^ 



(36) 



A constant value of d^Vl/dt^ implies that Vt{t) is a piece- 
wise quadratic polynomial. Whenever the second deriva- 
tive is reset, we choose integration constants such that 
VL and dfl/dt are continuous. Finally, we use n{t) to 
determine the actual rotation- matrix via Eq. (j34p . 

There are alternative control-feedback equations to 
Eq. (|36|) . Some of them are discussed in [40]. The de- 
tails of the feedback equation do not influence the main 
focus of this paper which is how to represent rotations 
and control parameters. 

In SpEC, Eq. ([M)) is integrated with a 5-th order 
Dormand-Prince time-stepper |48|. 

While Eq. (|34p analytically preserves the unit-norm 
of q, numerical integration will not identically preserve 
\q\ = 1. Therefore, the q{t) returned by the ODE- 
integrator is rescaled to unit- length, q —^ q/\q\ before 
it is used to construct rotations. 

Whenever d^n/dt^ is reset via Eq. dSS]), q of the ODE 
integrator is also rescaled to unit length. 



Equation (|35p can also be derived with the formal pro- 
cedure introduced in Section IIIBI This derivation will 
highlight an ambiguity not visible in Eq. (|35|) . and will 
also result in the control parameters for scaling and trans- 
lation. We start with 



X = aqxq* + T 



(37) 



where x,x,q,T are quaternions and a <E K. All vector 
quantities arc now treated as quaternions via the identifi- 
cation map V = (0, ii). We now perturb a — !■ a + Sa, T — > 

T+6T, q^qU^^^ 

in vectors v being mapped to 



^ 1 . The q - perturbation will result 



v' = q{l+^-S.yfl-^-S.)g*=q^g^ (38) 



where w = ll + -^)t)ll — j-j ■ This shows that the 

imaginary part Sq of 5q = (0, Sq) represents a rotation 
in grid coordinates. 

The quaternion version of Eq. Q is: 



q* +T + ST. (39) 



aqxiq* +T = 



with i ^ A,B. Because the real part of Eqs. (|39p are 
trivially satisfied, Eqs. p9p represent six equations, three 
each for black hole A and for black hole B. We seek to 
solve Eqs. ([Ml) for the unknowns da, ST ~ {0,5T), and 
Sq={0,dq). Because ST and Sq have three components 
each, we have in total seven unknowns. The additional 
degree of freedom arises because the rotation around C 
is not yet fixed. Recall that in the Euler angle repre- 
sentation, we remove this degree of freedom by setting 
4> ~ 0, c.f. Eq. ([7]). Expanding Eq. ([55]) to linear order in 
the perturbations and subtracting the equation for black 
hole B from that for black hole A, it is straightforward 
to show that 



Sa = 



X -C 



-la. 



- CxX ^ 

oq = — h ao. 



(40) 



(41) 



and 



(0, 5T) =aq[xA-CA-Sq^CA- —ca ] q* . (42) 



The normalization chosen corresponds to the fact that only the 
direction of the two vectors matter in the context of rotations, 
and due to the scaling control system we should have to first 
order, ||X|| ~ ||C||. 



In Eq. (U^l Sq A Ca = {0,Sq x ca) and Sq, Sa are to 
be substituted from Eq. (|i(Il HT|) . The parameter a in 
Eq. (PT|) is undetermined, reflecting the extra degree of 
freedom already mentioned after Eq. ([39]) . It parameter- 
izes the component of Sq parallel to C; i.e., a rotation 
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FIG. 2. Newtonian simulations with inclination of the or- 
bital plane of angle /3 = 0, 10, 70 degrees from the xy plane, 
performed with both control systems. Time is measured in 
units of orbital period. 




FIG. 3. Post-Newtonian simulations with inclination of the 
orbital plane angle /3 = 0, 10, 70 degrees from the x — y plane, 
performed with both control systems. Time is measured in 
units of initial orbital period. The binary is equal mass and 
non-spinning, with the initial coordinate separation of 20. 



about the axis C connecting the two excision spheres. 
We shall choose it to minimize the overall rotation \\5q\\: 



a = 0. 



(43) 



With this choice Eq. (HJ) simplifies to Eq. ^Q. The 
choice a = is equivalent to the minimal rotation frame 
of Boyle et al. [43] ; it minimizes artificial activity of the 
control system that is not connected to the physics of the 
binary black hole. 



III. NUMERICAL RESULTS 

To test our new approach to the rotation control sys- 
tem we begin with the simplest possible system that 
still exhibits the desired behaviour, namely a Newto- 
nian circular binary. We consider an equal mass, non- 
spinning, circular binary at separation of 20 M . The 
orbital plane is inclined with respect to the xy-plane by 
angles /3 = 0, 10, 70 degrees. The performance of the con- 
trol system is quantified by the magnitude of the control 
parameters \\Q\\ = \/ QiQ"^ where the summation extends 
over all the components of the control error for rotation, 
defined by Eqs. (|10b[ [TUcl) for Euler angles and Eq. (|iT|) 
for quaternions. Independent of the inclination /?, we al- 
ways initialize the control system as if the binary is in 
the xy-plane. This is of course only correct for /? = 0; 
for /3 ^ 0, the control system will also have to demon- 
strate that it can compensate for an utterly erroneous 
initialization. Figure [2] shows \Q\ for the 3 cases. For 
/3 = 0° both control systems perform very well with an 
extremely small value ofQ ~ 10~^^. For /3 ^ 0, there are 
initial transients due to the intentionally wrong initializa- 
tion of A^. These transients decay exponentially on the 
damping timescale of the control system; here, r = P/56 
where P is the orbital period. Once the transients have 
disappeared for /3 ~ 10°, the quaternion results are un- 
changed while the Euler-angle control error has increased 
by 6 orders of magnitude. Several periodic sharp features 
start to appear. Finally, for /3 = 70° the quaternion [|(3[| 
is again at 10~^^ while the Euler-angle ||Q[| grows by an- 
other two orders of magnitude and shows sharp oscilla- 
tory features, which ultimately makes the control system 
inviable. 

Figure [5] foreshadows already the main conclusion of 
this work: The Euler-angle approach depends on the 
plane of the orbit, and has increasing difficulty in con- 
trolling the coordinate mappings as the orbital plane be- 
comes orthogonal to the xy-plane. While the Euler-angle 
control system becomes singular only at exactly /3 = 90°, 
the effects of this singularity are already clearly visible 
for /3 = 10°. In contrast, the quaternion control sys- 
tem is rotationally invariant, and hence, it controls the 
coordinate mapping equally well for any inclination /3. 

Next, we turn to a more interesting test that also in- 
volves the control system for the expansion factor a(t). 
We consider a post-Newtonian equal mass, non-spinning 
black hole binary. The relevant FN equations of motion 
can be found in [50|. Figure |3] displays a set of three runs 
done with both control systems, again choosing to tilt 
the orbit relative to the xy-plane by angles /3 = 0, 10, 70 
degrees. 

After the initial transients due to intentionally wrong 
initialization of the control system, both of the con- 
trol systems handle the /3 = case equally well with 
1 1 (51 1 ~ 10~^ showing regular oscillations due to a smaU 
eccentricity of the orbit. When /3 = 10°, the quaternion 
control system performs in exactly the same fashion as 
for /3 = 0°. The Euler-angle system, on the other hand. 



d1 1 .68q2.5 




dll.68q2.5 



dll.68q2.5 



1000 



2000 3000 

Time (M) 



4000 



FIG. 4. The inclination angle, /3 for the three systems under 
study. 



begins to struggle: the amplitude of ||Q|| grows by about 
two orders of magnitude. The situation grows worse still 
for Euler angles when /3 = 70°, where the control error 
increases by another two orders of magnitude and sharp 
features appear. 

Meanwhile, the curves corresponding to the quaternion 
control system show exactly (to within numerical accu- 
racy) the same value of ||Q|| for all inclinations. This 
is exactly what we expect from a rotationally invariant 
control system - the orientation of the orbital plane is 
irrelevant. 

Finally, we test the quaternion-based control system 
using its main application: simulations of precessing bi- 
nary black hole systems in full numerical relativity. Quite 
generally, this precession may cause the orbital plane to 
rotate by 90 or more degrees with respect to the ini- 
tial conditions. The behaviour of the system depends on 
mass ratio and the two spin vectors. We choose a set of 
three simulations to be evolved using full numerical rela- 
tivity that exhibits mild to significant precession. Table 
|T] summarizes the initial conditions. The initial data was 
constructed from the superposition of two Kerr-Schild 
metrics for the conformal metric as in [51[ (so called SKS 
initial data). The eccentricity has been removed by an 
iterative process [S^], so that the final eccentricities for 
all three cases are a few xlO~^. 

Figure 2] shows the inclination angle f3 = 
arccos(f22/|f2|) which measures the angle between 
the normal to the instantaneous orbital plane and the 
initial direction of the normal, which is by convention in 
the z-direction. The high-frequency oscillatory features 
arc due to the nutation of the orbital angular momen- 
tum, while the secular evolution is due to precession. 
Notice that the dl4.5ql.5 run completes a full precession 
cycle at f = 400GM. The curves end when the black 
holes merge. The maximum inclination angles are 
similar to those used for the Post-Newtonian evolutions 
above. Figure [S] shows the trajectories of the black holes 
in inertial coordinates. 




FIG. 5. The trajectories of the centers of the apparent hori- 
zons of the black holes in inertial coordinates for the 3 simu- 
lations. Top to bottom: dll.68q2.5, dl2q2.5, dl4.5ql.5. The 
left panels show the projection onto the xy plane and the 
right, the xz plane 



Figure |6] presents \\Q\\ for the three runs done with 
the quaternion control system. In the main panel of Fig- 
ure [H] it is difficult to compare the ||Q|| of different runs 
because of the difference in orbital frequency. The in- 
set shows the same \\Q\\ residuals for the three simula- 
tions but timeshifted such that the orbital frequency of 
Ain = 0.025 occurs at t = 0. As one can see the con- 
trol error norms lie very close to each other, and exhibit 
qualitatively similar oscillations with virtually no sharp 
features. The remaining differences in behaviour are asso- 
ciated mostly with the different eccentricities as well as 
masses and spins. The growth of the control parameters 
with time is caused by the more rapid inspiral towards 
merger. Our numerical experiments demonstrate that 
the quaternion approach is indeed suitable for simulat- 
ing arbitrarily precessing configurations. 



Name 


q 


Xi 


X2 


Do/M 


do 


MQ.0 


dll.68q2.5 

dl2q2.5 

dl4.5ql.5 


2.5 
2.5 
1.5 


(0.000, 0.575, -0.556) 

(0.000, 0.410, -0.287) 

(0.000,0.285, 0.093) 


(0.000, 0.360,-0.347) 
(0, 0, 0) 
(0, 0, 0) 


11.68 

12 

14.5 


-0.000290589649 
-0.000108923113 
-0.000016947638 


0.02264246 
0.02180603 
0.01664958 



TABLE I. The initial conditions used for the numerical relativity runs. Given are the mass ratio q — mi/77i2, the dimensionless 
spin- vectors Xi ^nd Xa, the initial separation Do, initial radial velocity do, initial orbital frequency fio. The initial orbital 
angular momentum is in the z direction and the line connecting the two black holes is parallel to the a;-axis. 
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FIG. 6. Three full NR simulations performed with the quater- 
nion control systems. The initial conditions are listed in Table 

m 



IV. DISCUSSION 

Simulating precessing binaries poses a challenge for 
excision-based numerical techniques. This challenge is 
resolved here by developing coordinate mappings which 
make the black holes be at rest in grid coordinates. This 
transformation is dynamically controlled by a feedback 
control system since the trajectories of the black holes are 
not known in advance. In the most general case this map 
involves a rotation, and while Euler angle parametriza- 
tion works well for mildly precessing binaries, it exhibits 
coordinate singularities for polar orbits which leads to the 
breakdown of the simulation. To rectify the situation, we 
have created a control system that represents rotations 
using quaternions. Quaternions do not suffer from co- 
ordinate singularities and work for gcnerically precessing 
systems. The quaternion-based control system is able to 
successfully perform fully general relativistic simulations 
of highly-precessing binaries, allowing the investigation 
strongly precessing binary black holes and broadening 
the range of parameter space that can be explored. The 
techniques developed here have already been utilized in 
the simulations presented in [4l|, \^ 

The quaternion control system, as described and devel- 
oped above, is related to the minimal rotation frame [4^ 
(see also [5J]). For the control-system developed here, as 
in the minimal rotation frame, a preferred axis exists (the 
line connecting the two black holes vs. the instantaneous 



preferred emission axis of the gravitational waves). In 
both cases, the rotation about this axis is not a priori de- 
termined. And in both cases, this rotation is chosen such 
that the instantaneous rotation frequency of the rotating 
frame is minimal. In the present context, this condition 
is imposed by Eq. (^5]) . 

As a useful byproduct of the quaternion control sys- 
tem, one obtains an accurate estimate of the orbital fre- 
quency and the orbital phase during the numerical run 
without the need for any post-processing. The CI in 
Eq. p4p is the instantaneous rotation frequency of the 
grid frame relative to the inertial frame, given in compo- 
nents of the grid frame. Converting to the inertial frame, 



(0,(J) = qrjq* 



(44) 



If the control system were perfect - i.e. if Q = - 
then a; given by Eq. ([^^ would be the instantaneous 
orbital frequency. Because Q ^0, Eq. pi)) only gives an 
approximate orbital frequency, albeit a very good one: 
The upper panel of Fig. [7] shows the fractional difference 
between \lo\ from Eq. (|44[1 and the exact numerical orbital 
frequency obtained by post-processing. The difference 
oscillates around zero with relative amplitude of 1 x 10~'^. 
It is also straightforward to integrate 



= ir!| = 



(45) 



to obtain the orbital phase of the precessing binary. In 
practice, we add Eq. PSI) to the set of ordinary differ- 
ential equations Eq. (|34p that are integrated to obtain 
the rotation quaternion q{t). The difference between the 
orbital phase from the control system Eq. (P5|) and the 
exact orbital phase from the BH trajectories is shown in 
the lower panel of Fig.[7l The difference is ^ 10~^ radians 
until merger, during an inspiral lasting 105 radians. Inci- 
dentally, this again demonstrates that our control system 
works exactly as expected. 
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APPENDIX 

A. Rotation control parameters in matrix notation 

The underlying idea for the rotational control system 
that we have described in section lIlC3l is independent of 
the use of quaternions to represent rotations. For exam- 
ple, one could have used infinitesimal rotation matrices to 
achieve the same goal. Below is the demonstration of the 
same derivation as in section IIIC 31 but now in terms of 
a rotation matrix R. We start with the following version 
ofEq. (O: 

aRx^ + f^{a + 6a)RiI + 6R)c, + f + Sf (46) 

where as usual i = A, B, Xi,Ci,T, ST G M.^, I is the 
identity matrix, a, (5a £ R and R, SR G Max 3 . Once 
more we seek to solve this system of six equations for 
the unknowns Sa, ST, and SR. Note that since SR is an 
infinitesimal rotation matrix, it is skew symmetric and 
thus has three independent components, for a total of 
seven unknowns. Expanding Eq. (|46p to first order in 
perturbation and subtracting the equation for black hole 
B from that of black hole A one can show that: 
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SR, 



Sa 



<^ilkC 



X -C 



CxX 

\\cr 



aC, 



ST ~ aR \ XA — CA 



Sa ^ 

X CA CA 

a 



(47) 



(48) 



(49) 



These results match exactly Eqs. (gOll-dlll)- Thus we 
see that indeed infinitesimal rotation matrices could have 
been used to represent rotation. We selected quaternions 
for our work primarily for numerical reasons, the main 
being the ease of correcting numerical drift from a rota- 
tion, which for quaternions amounts to a simple renor- 
malization. 
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